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ABSTRACT 


This thesis describes the development of a Fortran IV 
computer program for finding shock-induced stresses in a 
Ere piping system. Discretized equations of 
motion are formed by the finite element method. Input 
motions are support translations specified by response 
shock spectra. Maximum responses in individual modes and 
resulting octahedral shearing stresses at selected points 
are found for shock motions in three orthogonal directions. 
Extreme stresses, based on the square root of the sum of 
the squares of the modal stresses, are estimated for each 
input direction. An example problem is analyzed to demon- 


strate the use of the program. 
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LIST OF SYMBOLS 


A single underline on a capital letter denotes a 


rectangular matrix, and a single underline on a lower case 


letter denotes a column vector. Superior dots denote time 


derivatives. The symbols used in the computer program are 


described in Appendix E. 
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Element cross-sectional area 

Mode participation factor, mode r 

A constant 

Pipe outside diameter 

Young's modulus 

Element generalized force vector 
Components of the generalized force vector 
Shear modulus 


Second moment of the pipe cross-sectional area 
about a diameter 


Identity matrix 

Mass moment of inertia per unit length about pipe axis 
Assembled stiffness matrix 

Element stiffness matrix 

Length of an element 

Assembled inertia matrix 

Element inertia matrix 

Element inertia submatrix, bending 1-2 plane 

Element inertia submatrix, bending 1-3 plane 
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Element inertia submatrix, torsional motion 
Modal mass, mode r 

Mass of the lagging per unit length 

Mass of the pipe per unit length 

Row vector of shape functions 

Vector of principal coordinates 

Erinclpal coordinate, mode r 


Moment of pipe cross-sectional area lying on one side 
of neutral axis 


Displacement vector 

Radius of gyration 

Base displacement 

Kinetic energy of an element 

Pipe wall thickness 

Vector of absolute displacements corresponding to 
unit base displacement 

Components of displacement vector 
Modal matrix 

Spectrum velocity, mode r 

Eigenvector, mode r 

Weight of lagging per unit length 
Weight of pipe. per unit length 

Absolute displacement vector 

One degree-of-freedom system coordinate 
Modified specific weight 

Dimensionless length coordinate 

Density 

Modal stress 


Shearing stress 


Octahedral shearing stress 


oct 

a? Spectral (diagonal) matrix 

ω Natural circular frequency, one degree-of-freedom 
system 

Wr Circular frequency of mode r 


Superscripts 


T Transpose of matrix 


(r) Mode designator for eigenvectors (r=1,2,...,n) 
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I. INTRODUCTION 


With the greater demand for electrical power, and the 
diminishing fossil fuel supply, there has been an increase 
faethe construction and шеге? nuclear power generating 
plants.  Correspondingly, there has been greater concern for 
the safety of the power plant during an earthquake. 

Miso, the reduction in the size of the Navy has ео 
an increased emphasis on the integrity of the internal 
systems of a naval vessel under shock loads. Thus, there 
has been a growing engineering interest in finding a rapid 
and accurate method of analysis to determine the shock- 
induced stresses in a complex continuous piping system. 

This thesis describes and presents a computer program, 
written in Fortran IV computer language, to find the 
Stresses in a piping system responding to shock. 

Complex piping systems can be analyzed by using dis- 
cretized models. In using the discretization technique, the 
continuous piping system is modeled by finite size elements 
connected at nodes. After the system has been subdivided, 
ΠΕ elastic and inertial matrices of each element can be 
found and assembled to form the elastic and inertial 
matrices of the system [1,2]. 

A considerable volume of material is available to assist 
in analyzing structures and systems responding to 
earthquakes [3,4,5]. Two of the techniques used in earth- 


quake response analysis are modal analysis with a 


ΠΠ 


discretized model and transfer functions [4,5]. The shock 
input Tor earthquakes can be Sspecii ted Ве ne 
history [3], or shock spectrum [4,5]. 

The Navy uses a technique called the Dynamic Design- 
Analysis Method (DDAM) to determine the shock response of 
shipboard equipment and systems [6]. The DDAM is a modal 
analysis technique using shock spectra to specify the shock 
inputs. 

The shock spectrum presents, as a function of the system 
natural frequency, the maximum response of a one degree-of- 
freedom system responding to a shock motion. The shock 
spectrum is generally presented in terms of response velocity 
ereacceleration. 

Several graduates of the Naval Postgraduate School have 
written theses on determining natural frequencies and mode 
shapes for piping systems. Fink [7] developed a program to 
analyze planar piping systems, including out-of-plane 
bending, using transfer matrices. Baird [8] presented the 
necessary theory to expand Fink's work to a three-dimensional 
piping system. Rudolf [9] developed a program that deter- 
mines the frequencies of a general three-dimensional piping 
system. Because the transfer matrix does not yield a dis- 
eretized model of the structure, there is no ready means 
for using the frequency and mode shape data to find shock- 
induced modal responses. On the other hand, the finite 


element technique provides a consistent discretized model 


за 


which may be used for finding frequencies, mode shapes, modal 

responses to shock inputs, and resulting stresses. For this 

reason the finite element method was chosen for this thesis. 
Shock spectra and mode participation factors are used 

to determine the modal responses. At selected points the 

modal octahedral shearing stresses are combined by scalar 

addition of the distortion energy to estimate extreme 


stresses. 
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TE, THEORY 


Material and structural linearity and material local 
homogeneity and isotropy are assumed in developing the 


mathematical model. 


A. MODEL LIMITATIONS 

There are few inherent restrictions on the capabilities 
of the finite element method to model complex details of 
piping configurations. Considerations of the quantity of 
input and output data, program length, core storage capac- 
ity, computing time, and cost do impose practical limita- 
tions. The limited time available in developing this 
thesis has necessitated the following restrictions on the 
model used here. 

The system consists of straight lengths of constant- 
section pipe (elements). Each element centerline is parallel 
to one of the three mutually orthogonal global refer- 
ence axes, and 90° bends are replaced by fictitious exten- 
sions of the tangent sections to the intersection point of 
the centerlines. Effects of added mass contributed by 
external lagging are represented, but no provision is made 
for added mass due to valves or fittings. Pipe hangers 
Furnishing uniaxial restraint in a global direction may 
also be represented. Bending action is described by the 


Euler-Bernoulli beam theory, neglecting shear deformation 


14 


and rotatory inertia. All ends of the coniieuration ane 
treated as fixed. The shock input motion is a translation 


of the base! along a global axis. 


В. FINITE ELEMENT METHOD 

Once the system has been subdivided into elements, 
moe clastic and inertial characteristics are determined. 
The element stiffness matrix can be determined by many 
different techniques, whereas the inertial matrix is gener- 
ally determined by using shape functions and integrating 
along the length, [1,2]. 

ΙΕ ΤΙ ΠΕΠ ΤΡΙ οεπσδας ο ποσο 

Consider an element with displacement components at 


the nodes (ends) as shown in Fig. 1. 





FIG. 1 


ELEMENT GENERALIZED DISPLACEMENTS 


the clamped ends of the configuration are treated as 
joined to a rigid base. 
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Here, the double arrowhead vectors are positive rotations 
in accordance with the right-hand rule. The displacement 


vector is 
“Т 
] 


prom Pig. 1, components ч, апа Ur are longitudinal dis- 


placements, components Чо, us» Up ue» Ug» бо» Uii and Чо 
are bending displacements and rotations in the 1-2 and 2-3 
planes, and components Uy and чуо are twisting rotations. 
Element Stiffness and Inertial Matrices 

Przemieniecki [2] forms the element stiffness ma- 
trix by solving the differential equations of the element 
displacements. Appendix A reproduces the element stiffness 
matrix with Przemieniecki's transverse shear deflection 
factor (¢) set equal to zero. 

To determine the element inertia matrix, an express 


sion for the element kinetic energy TŠ is formed, rotatory 


inertia due to bending being neglected. 


H 
11 
Mle 
Q 
>> 
о ~ 
| O 
= 
|= 
© 
Ω, 
ЧҮҮ 
11 
грн 
Kes 
|= 
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pere 4 is the element length, C is a constant, ΜΓ is the 
element inertia matrix, N is the shape function row Vectors 


and & is a dimensionless length coordinate. 
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Therefore 


ме = саг н! 


cs 


|= 


ЧЕ 


Since the local element axes 


(1) 


coincide with the cross 


section principal axes, the inertia matrix may be deter- 


mined from 2x2 and 4x4 element sub-matrices. 


Consider the 


following motions of the element: 


a) longitudinal: 
5 т 
ар, = lu, uy] 


апа С = 


element and p is the density. 


ме _ оА% e al 
—L 6 1022 


O) twist: 


2 T 
qp ^ [uy uj, Ny 


and C z 


length about the pipe axis. 


5 N 


= l — E, N. = € 


pA where A is the cross-sectional area of the 


From Eg. 1 


(la) 


= (1 ae ENR. N19 = oe 


J where J is the mass moment of inertia per unit 


From rg. 1 


(19) 


ο) bending in the 1-2 plane: 


_ T 
Чв12 = (42 96 98 912], М 


Nç = в(Е-2Е2+Е3), Ν 


== ге 3 2 
8 = 3Е -2t 5 N45 τ 


= 1-3&°+2Е&3, 


ее. 


n 


and C = pA. Prom Eg. 1 


15ο 7220995 -132 
= 2 2 
Maio ® 155 22% ar 130 232 (Te) 
54 132 156 -222 
-132 344 -22% ge 


Bending in plane 1-3 results in the same numerical coeffi- 
erents as MO ‚ Because the positive senses fon the 
rotation angles are reversed, rows 2 and 4 and columns 2 
and Y are multiplied by minus one. These submatrices are 
assembled to form the element inertia matrix. This matrix 


is displayed in Appendix A. 


C. ASSEMBLED EQUATIONS OF MOTION 
The governing equations of motion, in matrix notation, 


for undamped free vibration [2] with a fixed base are 


iS 


кже (2) 


where q is the assembled displacement vector and M and K 
are the assembled inertial and stiffness matrices. М апа 
pare referred to a global coordinate system and м® апа Ke 
ere referred to local coordinates. The assembly process 
consists of realigning the local system to the global 
System and building the M and K matrices. Reference [2] 


gives sample assembly techniques. 
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D. MODAL ANALYSIS 

Modal analysis is the process of finding the response 
motion of a many degree-of-freedom system by superposition 
of the response motions in the principal (or natural) modes. 
This technique is advantageous because the motions in the 
principal modes are independent of one another (uncoupled). 
To use the method, it is first necessary to find the nat- 
ural frequencies and mode shapes for the free motions 
governed by Eq. 2, i.e., to solve the eigenvalue problem 
K v = o^ M v. 

Finding eigenvalues ο) and eigenvectors (v) is a 
standard problem of numerical analysis. The first step is 
a transformation to a single-matrix problem. Details of 
the transformation used here, which requires a Cholesky 
decomposition of the inertia matrix, are given on p. 349 
of reference [1]. Following this, Jacobi rotations [15] 
are used to find the spectral and modal matrices (α΄ and V). 
The modal matrix is normalized with respect to the inertia 


T 


КОШ Лу, Е, М ш у = Lo where, Ll is, the nth order identity 


matrix. 

The shock analysis is based on using velocity shock 
spectra to characterize the input motion and mode parti- 
cipation factors to determine the responses. Although both 


shock spectra and mode participation factors have been 


extensively used in earlier studies [6], the finite element 
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discretization necessitates a new development of working 
formulas. Details of this development are given in 
Appendix B. 

Each shock input motion consists of base translation 
parallel to a global axis. It is shown in Appendix B that 
the mode participation factor 5, for mode ο rs 


ως 


b =v Mu | (8.5) 


5 — 


(т) 


лете у is the modal eigenvector and u is спе тес оте 
absolute displacements corresponding to unit base displace- 
ment. The maximum relative displacements due to mode r 
response are 
qo vU). VL (В. ба) 
—max — IE CAD 
where Wn äsr che medalzerreurarzıı zequeneyzaınd V. is the 
Spectrum velocity at that frequency. 

Since responses to shock inputs in each of the ο. 
epal directions are Separately determined, the vector u 
and the mode participation factors br must be evaluated 
separately for each direction. The modal masses m, (see 


Appendix B) are likewise different for each input direction. 


E. GENERALIZED FORCES AND STRESSES 
The subvector 5 of element peak modal displacements 


may be found from the system vector, taking into accouny cmc 
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orientation of the local reference axes. From this the 
element generalized force vector τι is determined. This 
force vector consists of the elastic and inertial 


contributions for each element. 
e e e 
f = (К-ы М) а ο) 


Fig. 2 shows the positive directions for the components 
of the generalized force vector. 
71 


козе га Sm, {12 





πα. 5 


ELEMENT GENERALIZED FORCES 
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The stresses at the nodal points of the element can be 
determined from the generalized forces by standard methods 


ӨГ solid meehantcs [19]. 


ΕΞ 





ТТТ. PROGRAM DEVELOPMENT 


The computer program was written in Fortran IV com- 
puter language using double-precision arithmetic on an 
IBM 360/67 digital computer. It consists of a main calling 
program and seven Subroutines. The main program calls the 
eu routines, where the actuari caleulations occur, in the 
necessary order. Information is passed between the main 
program and the subroutines by the use of a common storage 


соге. Appendix E contains the program listing. 


A. SYSTEM DESCRIPTION 

After the piping system has been modeled and subdi- 
vided, subroutine INPUT is used to read the data input of 
the geometry and material properties of the system. The 
Beemetry of the piping system model is determined by the 
global coordinates of the nodes of the model. The global 
system is a right-handed orthogonal coordinate system. 
Bach pipe segment must be parallel to a global axis, but 
there are no restrictions on the placement of the origin of 
the global system. The parameters from the subdivision of 
the system are read. For each node the global coordinates 
and boundary conditions are read. If the node is fixed, 
the stiffness and inertial matrix components due to the 
node displacements are not assembled into the system 


matrices. 
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The numbering of nodes and elements is arbitrary. 
There are two nodes per element and six degrees of freedom 
per node. For each element, the node numbers and pipe 
group are read. A pipe group includes all elements that 
have the same Young's modulus, Polsson's ratio, specific 
weight, radius of gyration, and cross-sectional dimensions. 
Appendix D shows the method used to calculate the effect 
of lagging on the specific weight and the radius of 


gyration. 


B. STIFFNESS AND INERTIAL MATRIX FORMULATION 

Subroutine FØRM generates the element stiffness and 
inertia matrices as shown in Appendix A. Rotatory inertia 
and shear deformation due to bending are neglected. The 
orientation of the element is determined with respect to 
the global coordinate system and the correspondence between 
local and global degrees of freedom is established. After 
determining the correspondence between the local and global 
degrees of freedom, the element stiffness and inertia 
matrices are assembled to form the system stiffness and 
lnertia matrices. If pipe hangers are present, the corres- 
ponding (uniaxial) hanger stiffnesses are added to the 
appropriate diagonal elements of the system stiffness 


matrix. 


C. FREE VIBRATION MODE SHAPES AND FREQUENCIES 
Subroutine CHÓMÓD uses Cholesky decomposition of the 


system inertia matrix and coordinate transformation to form 
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a single symmetric matrix for which eigenvalues and eigen- 
vectors are found. Subroutine JACRØT is a modification of 
a program originally coded by Professors G Cantin. το 
uses the Jacobi variable threshold method to find eigen- 
values and eigenvectors of a real symmetric matrix. The 
modal matrix, when transformed back to system coordinates, 
is normalized with respect to the system inertia matrix. 
Since the system matrices increase in size with finer 
subdivisions or more complex systems, an economizing tech- 
nique, to reduce the number of degrees of freedom in the 
eigenvalue problem, was investigated during the program 
development. The component mode synthesis method as used 
by Benfield and Hruda [11] was studied. It was established 
that the constrained component branch technique [11], with 
interface loading, with no node suppression gave results 
identical with those obtained by the direct method described 
above when applied to planar vibration of a clamped-clamped 
beam represented by six elements. Despite the success of 
this trial, it was concluded that the added program com- 
plexity accompanying the use of component mode synthesis 
would offset the potential economies. Another standard 
economizer technique [14] was considered, but ultimately 


rejected for the same reason. 


D. MODAL RESPONSE 
The number of modes to be used in the stress analysis 


and the specification of the shock spectrum are input data 
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for subroutine SPCTRM. There are three choices available 
to select the number of modes to be used: a designated 
oper frequency limit, an upper limit whiten 13 а пе те 
of the fundamental frequency, and an option for any method 
Fie user desires. Three choices are also available for 
specifying the shock spectrum: a spectrum defined by 
straight-line segments; a constant velocity, then constant 
acceleration shock spectrum; and an option for any method 
ame user desires.  Thewehock input 15 restricted to a 
translation of the base along a global axis and is speci- 
fied by one shock spectrum and three scaling factors for 
the three input directions. Subroutine MØDE modifies the 
spectrum velocity for modal mass, if desired. There are 
three choices available: no correction, log of the modal 
Mergent correction, and an option for any method the user 
desires. This subroutine also calculates the mode 


mercicipation factors. 


E. STRESSES 

Subroutine STRESS calculates stresses at the nodes of 
the elements. The element peak displacements are deter- 
mined and realigned from global degrees of freedom back To 
Ehe local element degrees of freedom. The element gener- 
alized force vector is then caleulated. Consider the 
action of the generalized force vector on the element 


snown in Fig. 3. 
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υπο > 
STRESS LOCATION POINTS 


locations a,b,c, and d, the stress vector acting on the 
emess-section may be resolved into normal and shearing 
stresses. Fig. 4 shows the positive sign convention used 
for the normal and shearing stresses. 

The normal stress o and the shearing stress ı at each 
point shown in Fig. 3 are given by: 


1. at location a) 


0/21 


а 
I 


-f,/A - f. 


A 
I 


al 


t 


FIG. 4 
POSITIVE STRESS CONVENTION 


un re I is the second moment of cross-sectional area about a 
diameter, D is the outside diameter, Q is the first mo- 

Menu, about a diameter, of the portion of the cross-sectional 
au iying on one side of the diameter, and t 15 the pipe 
wall thickness. 


2 at location o) 


а 
tl 


£,/A + fii D/21 


a 
II 


fig D/AI - fg 0/2 16 


ay at location Б) 


o = - f1/A - fe DOT 


Q/2 It 


A 
I 


ἘΠΕ D/AI + f, 
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Ш ас location d) 


UM NE Ε/Α - - D/2T 


τ D/4I - 8, Q/2It 


£10 9 


Mise normal and shearing stresses are then used to find 
the octahedral shearing stresses [12], which are given 
by 


1 
Т τὸ σος (4) 


oct 
For each element the above calculations are performed for 
each mode. At each of the four points, the shearing 
stresses for the individual modes are then combined by 
determining the square root of the sum of their squares. 
ES process is accomplished for each input direction. 
mis is done for each element in turn, starting with the 


"ρου. 


ГОО PROGRAM OUTPUT 

All input information is echo-checked. The square 
of the mode frequency and the mode shape are printed for 
each mode. The spectrum velocities and mode participation 
factors are also printed. The octahedral shearing stresses 
at points a, b, c, and d (Fig. 3) of each element are 


printed for the three shock input directions. 
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IV. PROGRAM TESTING 


In order to explore the capabilities of the 
and verify its integrity, a number of plane test 
ations were studied. These are described in the 


sections. 


fe (TEST PROBLEM 1 
The configuration shown in Fig. 5 has three 
uniform runs of equal length between the central 


and the clamped edges. 


πια Б 
TEST CERRO SN EDEN i 


program 
confipur- 


following 


identical 


Junetdionm 


The mode shapes and frequencies for in-plane vibration of 


this system were studied extensively with a developmental 


program. The present program gave identical results for the 
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in-plane modes. The out-of-plane modes exhibited the 
required symmetric or anti-symmetric form. 

Additional program tests on this configuratlon involved 
Gl ferent element nuúnmbening and node numbering andic πα 
different orientations of the global axes relative to the 
structure. Mode shapes and frequencies were unaffected by 


these changes. 


Be) TEST PROBLEM 2 
This system, shown in Fig. 6, consisted of a straight 
ПП ουπ pipe clamped at the ends and represented by four 


elements. 


&— f 


FIG. 6 


TEOT PIPING SYSTEM 2 


The frequencies found for the lower modes showed good 
agreement with exact results. Table 1 shows the comparison 
of the first three bending frequencies of the finite element 
solution with exact results [13]. All the eigenvectors 
showed the required symmetric or anti-symmetric form. 
Bending modes occurred in pairs having equal eigenvalues 

and with the corresponding eigenvectors representing 


eerleetrions in Orthogonal planes. 
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a found for this system were studied in detail. 
At each node, the stresses are calculated at two poimes on 
the cross-section. Except at the clamped ends, there are 
two adjacent elements sharing each node, so that stresses 
at these sections are calculated twice. Complete agreement 
was found between these two sets of stress evaluations. 
miso, identical stresses resulted from shocks in the two 


directions perpendicular to the length. 


TABLE I 


COMPARISON OF FIRST THREE BENDING FREQUENCIES OF TEST 
PIPING SYSTEM 2 WITH EXACT RESULTS 


Mode Exact Finite Percent 
Element Difference 
(Hz) (Hz) 
1 422.32 422.65 0.08 
2 1161715 1171226 0.9 
3 2282.20 2329.64 cT 


C. TEST PROBLEM 3 

The equal-legged configuration of Fig. 7 consists of 
uniform pipe throughout. Mode shapes of this system again 
showeä the required symmetric or anti-symmetric form. 
Likewise, the shock-induced stresses exhibited the expected 


Symmetry properties. 
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Вто. Г 
TEST ΕΕ УЬТЕЙ 3 
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У. EXAMPLE PROBLEM 


ДО PIPING SYSTEM 

The example system analyzed is a model of the main 
steam piping system on a modern naval vessel. The main 
steam piping from the boiler to the rigid cross-connect 
anchor assembly is modeled. The ends of the piping system 
are clamped, and there are two hangers. The system is 


three-dimensional with no branches. The pipe properties 


are: 
outside diameter 7.625 inches 
inside diameter 6.011 inches 
Poisson's ratio 0.300 
Young's modulus 30x10 psi 
specific weight 0.327 1b/in? 
radius of gyration 4.04 inches 


The specific weight is a fictitious value which accounts 
for the weight of five inches of lagging bonded to the pipe. 
ВЗР. 8 15 а schematic of the piping system. The piping 

section Pl is 84 inches long, P2 is 264 inches long, P3 is 


108 inches long, and P4 is 288 inches long. 


В. DATA INPUT AND OUTPUT 

The input data cards are punched in accordance with the 
instructions contained in Appendix E. The echo-check of 
the data is shown in Table II and is arranged as it would 


appear at the end of the program deck. The system 
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Ей. 8 


EXAMPLE PIPING SYSTEM 


eigenvalues and eigenvectors are listed in Appendix C, 
Table VIII. Tables III and IV show the modal velocities 
aná the mode participation factors. The stress output 15 
given in Table V. 

This problem required a storage capacity of 174K bytes, 


and 2.67 minutes of computer time. 
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DEG OF FREEDOM 
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NUMBER OF PROBLEMS TO 8E SOLVED 
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TABLE II DATA INPUT, EXAMPLE PROBLEM 
ELEMENTS NODES 
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NODE 
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ELOCITY 
100.00 
1.000 


1 
Е V 
NCHES/SEC 


FREQUENCY 
1.000 _ 
MODE WEIGHT CORRECTION TYPE 


59.00. 
SHOCK INPUT FACTORS 


1.000 





TABLE 111 MODE VELOCITIES EXAMPLE, PROBLEM 
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TABLE V OCTAHEDRAL SHEARING STRESSES, EXAMPLE PROBLEM, PSI 
| ELEMENT 


ELEMENT 
9.27155D 
6.67426D 
6.60959D 
4.58857D 


ELEMENT 
6.60959D 

` 4588570 
1.879580 
4.960030 


ELEMENT 
5.164590 
1.888130 
4.458580 
1.114120 


ELEMENT 
4.45858D 
1.114120 
4. 561880 
7.499480 


EL EMENT 
4.561880 
7.499480 
4.543860 
6.252800 


ЕСЕМЕМТ 
4.543860 
6.252800 
3.613900 
1.086730 


1 
03 
03 
03 
03 


03 
03 
04 
03 


03 
04 
03 
04 


03 
04 
03 
03 


03 
03 
03 
03 


03 
03 
03 
04 


6.508980 
1.687530 
6.198620 
8.645530 


6.198620 
8.645530 
T.52TTTO 
7.392360 


7.613430 


6. 471280 


5.761720 
4. 914690 


5. 761720 
4.914680 
9.146560 
4. 300960 


9. 146560 


4.300960 


7. 596550 
4.330590 


7.596560 
4. 330590 
4.921410 
4.381410 


03 
04 
93 
03 


03 
03 
03 
03 


03 
03 
03 
03 


93 
03 
03 
93 


03 
03 
03 
03 


03 
03 
03 
03 


8.408350 
6.482560 
2.617210 
3.095510 


2.617210 
3.095510 
7.714060 
1.048080 


1. 776200 
7.523890 
3.033880 
2.613220 


3.033880 
2.613220 
3.568610 
6.189240 


3.568610 
6.189240 
4^. 165240 
7.081600 


4.765240 
7.081600 
2.472250 
7. #11130 


03 | 
03, 
03 
03. 


Done АМА: а асаа е =. 


1.055150 
4.203340 
5.886870 
4.312620 


ELEMENT 
5. 886870 
4.312620 
2.656300 
7.424390 


© ELEMENT 


03 
03 
03 
03 


03 
03 
03 
03 


03 
03 


03. 


03 


03 
03 
03 
03 


7.592110 
2.46920D 
5. 176440 
4.857170 


ELEMENT 
5.716440 
4.857170 
7.218960 
9.868170 


ELEMENT 
6.940940 
9.868170 
1.043000 
3.535350 


. ЕСЕМЕМТ 


+ 


1. 943020 
3.535340 
1.769770 
6.153500 


- = rep 
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T 
04 
03 
03 
03 


93 
03 
03 
03 


03 
03 
93 
03 


10 
03 
03 
03 
03 


11 
03 
93 
04 
93 


12 
94 
03 
04 
03 


3.195620 
5.395320 
3.140880 
5.552520 


3. 140880 
5.550520 
3.213450 
1.318900 


1.272420 
2.811240 
1.912420 
2.063100 


1.010420 
2. 063100 
6.128630 
3.488530 


6. 083360 
3.488530 
4.319220 
2. 305370 


4. 319220 
2.305370 
1.275160 
2. 751500 


03 
93 


03 


03 


эз 
03 
03 
94 


04 
03 
04 
03 


04 
03 
03 
03 


03 
93 
03 
03 


33 
03 
04 
03 


1.845000 
2.892360 
3.869280 
3.820040 


3.869280 


3.820040 
3. 5677150 
4.405200 


4.365370 
3.508090 
3.201110 
8.551300 


3.201110 
8.551300 
4.520530 
1.570740 


4.585790 
1.579740 
2.335150 
9.553940 


2.335750 
9.553940 
6. 004 TTD 
1.630620 


03 
03 
03 
03 


93 
03 
03 
03 


03 
03 
93 
03 


03 
03 
03 
04 


03 
04 
03 
03 


03 
03 
03 
04 





Vi.» DISCUSSION 


Some of the limitations of the capabilities of the 


program developed above are considered here. 


A. CONFIGURATION LIMITATIONS 

The restrictions that pipe axes must be parallel to a 
global axís and that finite radius bends are not represented 
provide the principal limitations on the configurations that 
δε modeled. There are no inherent limitations on allow- 
able topological complexity. 

Certain additional features of real piping systems that 
are not represented in the present modeling include added 
feos due to fittings and pipe contents, and partial fixity 
or elastic restraint at ends or intermediate points. 

Despite these limitations, it is believed that a signi- 
ΙΙ. fraction of current piping systems can be adequately 


modeled using the present program. 


B. SIZE LIMITATIONS 

For the present purpose, the appropriate measure of 
οτε size is the number of degrees of freedom of the 
system model. This number n, which bears no direct relation 
to physical size, determines the core storage requirements 
and execution time of the program. For approximate estima- 
tion, storage requirements are proportional to ae and 

5 


execution time is proportional το п”. Using the data from 


the example problem, one can estimate that a 100 
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degree-of-freedom system would require about HOOK bytes of 
core storage and about 9 minutes execution time. 

Increasing the problem size also increases the round-off 
errors. Because double-precision arithmetic (56 bit mantissa) 
is used throughout, observed round-off effects have been 
found negligible (a maximum of 1 unit in the 6th significant 
digit for n = 66) in all applications to date. The very 
small errors detected are attributed to the residual eigen- 
vector errors present "m the Jacobi rotations are 
terminated. 

In view of the foregoing considerations, it is believed 
that the practical upper limit on problem size (with the 
IBM 360/67) is about n = 110 and is determined by core 


storage capacity. 
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VII. CONCLUSIONS 


It is concluded that an effective program has been 


developed for determining shock-induced stresses in piping 


systems. The principal limitations of the present version 


can be removed by adding features that are clearly within 


the current state-of-the-art. Recommended additions are 


listed below: 


Hh 


2. 


Remove the pipe axes orientation restriction so that 
a general piping system can be analyzed. 


Develop. element stiffness and inertia matrices for 
bends. 


Modify the program so that partial fixity and elastic 
restraint at the ends or intermediate points may be 
included in the boundary conditions. 


Include the mass effects of valves, fittings, and 
pipe contents. 


Replace the Jacobi rotation method by a more efficient 
elgenvalue-eigenvector algorithm. 
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АРРЕМОТХ А 
ELEMENT STIFFNESS AND INERTIA MATRICES 


TABLE VI Element Stiffness Matrix 
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o о -® o ο ο ο ο ΞΕ 
= 7 τ Т 
оши о оф A Ne = 
BE E "p Ў 
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APPENDIX B: VELOCITY SHOCK SPECTRA AND MODAL RESPONSE” 


A shock spectrum exhibits the maximum response dis- 
placement of a single degree-of-freedom system whose base 
is subjected to the shock motion. Consider & single degree- 
of-freedom whose base has a shock displacement S, and whose 
displacement relative to the base is Z. The equation of 
motion is 


z i = - 5 (В.1) 


Where 15 the natural circular frequency of free vibration 
Deich the base fixed (s = 0). “If глад represents the ex- 
treme value of z in response to the shock motion s, then the 
displacement shock spectrum for that shock motion is a plot 
OT ας versus w. Equivalent information can be given in 
Mac storm of a velocity shock spectrum. In this form, the 
spectrum velocity V is related to the maximum displacement 


2 max by 


<: 
l! 


Q sas (Bez) 


То develop the equations for finding modal response 


"πο shock spectrum Gata, ΤΕΝ δε πε vector Of absolute 


lmhis Appendix is based on material presented by Prof. 


R.E. Newton in the course ME4522, September 1971. 
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displacements in the system degrees-of-freedom. For a 
uniaxial translation s of the base, this may be expressed 


as 


Ww coge us (B. 3) 


where q is the vector of displacements relative to the 
Base, and u is a vector of constants. To illustrate the 
meanings of these vectors, consider the beam element of 


EXE. 9. 





NNN Inertial Reference SEN 


FIG. 9 


BEAM ELEMENT COORDINATES 


À P T 
For this example, q = Ca, αρ 43 а] . ΝΞ [91 + 9, q4*s m > 


and u = [1 0 1 0]. It can be seen that u represents the 
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absolute displacements resulting from unit base 
translation. 
The equations of motion of an n degree-of-freedom 


system with base motion may be written 


W+XKg=0 


[3 


Using Eq. B.3, this becomes 


а+Ка = - Миз 


|= 


ШЕ the substitution q = V p where p is the vector of 


principal coordinates, this may be rewritten as 


VMVptVEKVp--VMus 
ΟΥ 

.. 2 .. 

Е O Ὁ 5 (B.4) 
шеге b = йт Mu. Eq. В.Д is equivalent to the scalar 
equations 

p.*u^p --b s, (r=1,2 n) (B. lla) 

T P τ D А нь : 


where Pr and ο are the rth components ог р rand be 


respectively, and Wn is the circular frequency of mode r. 
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The component ры is called the mode participation factor 
вот mode r. it is “given by 


T 
b. = yT) Mu (B.5) 


r 
where ντ) is the eigenvector for mode r (rth column of 
v). 

Comparing Eq. B.4a with Eq. B.1, and using Eq. B.2, 


the maximum response in mode r may be expressed as 


(po) sas b u Velas (B.6) 
where M is the spectrum velocity at frequency ш. Using 
Eq. B.6, the corresponding extremum for the relative 
displacement vector q is 

CEU | 

Qe ov bo Tu. (В.ба) 
Eq. B.6a is used for calculating peak modal responses ШЙ 
order to find stresses. 

In shipboard shock studies, it has been found that 
massive equipment may partially suppress the base motion 
[6]. This may be taken into account by appropriate modifi- 
cation of the spectrum velocity u Such ἃ correction τ᾽ 
quires apportionment of the system mass among the modes. 

The modal mass for the rth mode is taken to be 
2 


nes EA (B.7) 
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APPENDIX C: EXAMPLE PROBLEM EIGENVALUES AND EIGENVECTORS 


Table VIII shows the first six eigenvalues and the 
corresponding eigenvectors of the example problem. The 
eigenvalue (square of the circular frequency) appears at 
the head of a column followed by the components of the 
eigenvector grouped by nodes. The remaining eigenvalues 
and eigenvectors are omitted. Spacing of the eigenvalues 
remains approximately uniform throughout the remainder of 


Ber. The largest elgenvalue 3s 1.32 x 109. 
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TABLE VIII 


EXAMPLE PROBLEM EIGENVALUES AND EIGENVECTORS 
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APPENDIX D: CALCULATION OF MODIFIED SPECIFIC WEIGHT 


AND RADIUS OF GYRATION 


There is negligible effect on the stiffness of the 
pipe element from the lagging on the pipe; however, the 
lagging does contribute significantly to the mass of the 
element. By appropriately modifying the pipe element 
specific weight and radius of gyration, the mass effect of 
the lagging can be included. The modified specific weight 
is 


puc 2 (0.1) 


where De LS VG modified specific weight of the pipe ele- 
ment, and the sum us = Ντ) is the combined weight of the 
lagging and pipe per unit length. The radius of gyration 
is given by the relation 


2 


> το = (m + тү) n (6229 


L p 
where Јр апа Jr are the mass moment of inertia per unit 
length of the pipe and lagging, respectively, (m, + тү) 15 
the combined mass of the pipe and lagging per unit length, 


wda r is the radius or pyratiomn. 
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